# Script for Mallinson, Daniel J. 2014. "Upstream influence: The positive impact of PAC contributions 
#				on Marcellus Shale roll call votes in Pennsylvania"
#		  		Interest Groups & Advocacy


rm(list = ls(all = TRUE))

library(foreign)

#Read dataset for analysis
marcellus <- read.dta("upstream_data.dta")
data <- marcellus
data <- na.omit(data)

#################################  Descriptives for Appendix A ########################################

descriptives <- data[c("donate0810", "republican", "ideology", "unemploy", "county_rev", "margin", "senate", "tenure_total")]
sapply(descriptives, mean, na.rm=TRUE)
sapply(descriptives, sd, na.rm=TRUE)

################################ Analysis for Appendix B  #############################################


#Model of all votes
install.packages(c("robust", "sandwich"))
library(robust)
library(sandwich)
model.all <- glm(industry_vote ~  donate0810 + republican + ideology + unemploy + county_rev + margin + senate + tenure_total, family=binomial(link="logit"), data=data)
cov.all <- vcovHC(model.all, type="HC0")#Caluclate robust standard errors
se.all <- sqrt(diag(cov.all)) 
r.est.all <- cbind(Estimate = coef(model.all), `Robust SE` = se.all, `Pr(>|z|)` = 2*pnorm(abs(coef(model.all)/se.all), lower.tail=FALSE), LL = coef(model.all) - 1.96*se.all, UL = coef(model.all) + 1.96*se.all)
r.est.all

#Model of final votes
final <- subset(data, final==1)
model.final <- glm(industry_vote ~  donate0810 + republican + ideology + unemploy + county_rev + margin + senate + tenure_total, family=binomial(link="logit"), data=final)
cov.final <- vcovHC(model.final, type="HC0") #Caluclate robust standard errors
se.final <- sqrt(diag(cov.final))
r.est.final <- cbind(Estimate = coef(model.final), `Robust SE` = se.final, `Pr(>|z|)` = 2*pnorm(abs(coef(model.final)/se.final), lower.tail=FALSE), LL = coef(model.final) - 1.96*se.final, UL = coef(model.final) + 1.96*se.final)
r.est.final

#Model of amendment votes
amendments <- subset(data, amendment==1)
model.amend <- glm(industry_vote ~ donate0810 + republican + ideology + unemploy + county_rev + margin + senate + tenure_total, family=binomial(link="logit"), data=amendments)
cov.amend <- vcovHC(model.amend, type="HC0") #Caluclate robust standard errors
se.amend <- sqrt(diag(cov.amend))
r.est.amend <- cbind(Estimate = coef(model.amend), `Robust SE` = se.amend, `Pr(>|z|)` = 2*pnorm(abs(coef(model.amend)/se.amend), lower.tail=FALSE), LL = coef(model.amend) - 1.96*se.amend, UL = coef(model.amend) + 1.96*se.amend)
r.est.amend

#Make LATEX tables
install.packages("apsrtable")
library(apsrtable)
apsrtable(model.all, model.final, model.amend , se=c("vcov"), digits=3, stars="default")

#Run Wald test to calcualte chi-sq for comparing to original model
install.packages("aod")
library(aod)
wald.test(b=coef(model.all), Sigma=vcov(model.all), Terms=2:9)
wald.test(b=coef(model.final), Sigma=vcov(model.final), Terms=2:9)
wald.test(b=coef(model.amend), Sigma=vcov(model.all), Terms=2:9)



#############################################################################
########################## Figures ##########################################
#############################################################################

########################### Figure 1  #######################################
annual <- read.dta("annual_donations.dta")
time <- ts(annual$total, start=2000, frequency=1)

#jpeg("figure1.jpg", width=6, height=6, units="in", res=800)
pdf("figure1.pdf", width=6, height=6)
xrange <- c(1999,2012)
yrange <- c(250000,1700000)
ticks <- pretty(yrange)
plot(annual$year, annual$total, ylab="Total Annual Donations", xlab="Year", xlim=xrange, ylim=yrange, xaxt="n", yaxt="n", pch=20)
lab1 <- c("$502,802", "", "$535,003", "", "$842,271", "", "$810,718", "", "", "", "$1,593,372", "")
lab2 <- c("", "$290,944", "", "$278,269", "", "$301,258", "", "$456,506", "", "", "", "$560,861")
lab3 <- c("", "", "", "", "", "", "", "", "$540,876", "", "", "")
lab4 <- c("", "", "", "", "", "", "", "", "", "$1,091,215", "", "")
text(annual$year, annual$total, labels=lab1, pos=3, cex=0.7)
text(annual$year, annual$total, labels=lab2, pos=1, cex=0.7)
text(annual$year, annual$total, labels=lab3, pos=4, cex=0.7)
text(annual$year, annual$total, labels=lab4, pos=2, cex=0.7)       
lines(annual$year, annual$total)
axis(1, at=c(2000:2012), labels=TRUE)
axis(2, at=ticks, labels=c("$0", "$500,000", "$1,000,000", "$1,500,000", "$2,000,000"))
dev.off()


########################### Figure 2a #######################################

legis <- read.dta("house_senate_donations.dta")
over0 <- subset(legis, total>0)

#t.test for donations by party
t.test(over0$total~over0$republican)
#t.test for donations by shale county
t.test(over0$total~over0$shale)

#Crosstab of donations by party and shale district           
table.donations <- xtabs(over0$total ~ over0$republican + over0$shale)
summary(table.donations)
prop.table(table.donations)
#Crosstab of districts by party and shale district
table.districts <- xtabs(~ over0$republican + over0$shale)
summary(table.districts)
prop.table(table.districts)

#Mosaic of districts by party and shale county               
rep.shale <- 62
dem.shale <- 41
rep.no <- 52
dem.no <- 46
prop.rep.shale <- sprintf("%1.0f%%", 31)
prop.dem.shale <- sprintf("%1.0f%%", 20)
prop.rep.no <- sprintf("%1.0f%%", 26)
prop.dem.no <- sprintf("%1.0f%%", 23)
                 
prop.vec.mosaic <- c(dem.shale, rep.shale, dem.no, rep.no)
prop.matrix.mosaic <- matrix(prop.vec.mosaic, nrow=2, ncol=2, byrow=T, dimnames=list(c("Yes", "No"), c("Democrat", "Republican")))
prop.table.mosaic <- as.table(prop.matrix.mosaic)

#jpeg("figure2a.jpg", width=6, height=6, units="in", res=800)
pdf("figure2a.pdf", width=6, height=6)
par(mar=c(.5,3,2.7,1.2))
mosaicplot(prop.table.mosaic, color = c("gray60", "gray80"), main = "", shade = FALSE, xlab = "", ylab = "", las = 1, cex.axis =.8)
text(.36,.28,"62 Districts") #add counts by hand (note: need to use trial and error for coordinates)
text(.36,.22, "(31%)")
text(.36,.80,"41 Districts")
text(.36,.74, "(20%)")
text(.78,.28,"52 Districts")
text(.78,.22, "(26%)")
text(.78,.80,"46 Districts")
text(.78,.74, "(23%)")
mtext("Legislator Party", 2, font = 2, line = 1, cex = 1.1)
mtext("Presence of Shale Gas", 3, font = 2, line = 1, cex = 1.1)
dev.off()

########################### Figure 2b #######################################

rdonate.shale <- 68
ddonate.shale <- 19
rdonate.no <- 11
ddonate.no <- 2
                 
donate.prop.vec.mosaic <- c(ddonate.shale, rdonate.shale, ddonate.no, rdonate.no)
donate.prop.matrix.mosaic <- matrix(donate.prop.vec.mosaic, nrow=2, ncol=2, byrow=T, dimnames=list(c("Yes", "No"), c("Democrat", "Republican")))
donate.prop.table.mosaic <- as.table(donate.prop.matrix.mosaic)

#jpeg("figure2b.jpg", width=6, height=6, units="in", res=800)
pdf("figure2b.pdf", width=6, height=6)
par(mar=c(.5,3,2.7,1.2))
mosaicplot(donate.prop.table.mosaic, color = c("gray60", "gray80"), main = "", shade = FALSE, xlab = "", ylab = "", las = 1, cex.axis =.8)
text(.52,.37, "68%")
text(.52,.88, "20%")
text(.94,.37, "11%")
text(.94,.90, "2%")
mtext("Legislator Party", 2, font = 2, line = 1, cex = 1.1)
mtext("Presence of Shale Gas", 3, font = 2, line = 1, cex = 1.1)
dev.off()


########################### Figure 3a #######################################
var.name <- c("intercept", "Donations","Republican", "Conservative", "Unemployment", "Revenue", "Margin", "Senate", "Tenure")
or <- rbind(cbind(model=1, OR = exp(r.est.all[,1]), lci = exp(r.est.all[,4]), uci=exp(r.est.all[,5])), cbind(model=2, OR = exp(r.est.final[,1]), lci = exp(r.est.final[,4]), uci=exp(r.est.final[,5])), cbind(model=3, OR = exp(r.est.amend[,1]), lci = exp(r.est.amend[,4]), uci=exp(r.est.amend[,5])))

var <- c(var.name, var.name, var.name)
or <- data.frame(var, or)

or <- or[order(or[,1], or[,2]),]

#Delete intercept
or <- or[-7,]
or <- or[-7,]
or <- or[-7,]

nr <- nrow(or)

#Plot OR for all other covariates
or.norep <- subset(or, var!="Republican")
or.norep <- subset(or.norep, var!="Senate")
or.norepplusna <- lapply(split(or.norep, or.norep$var, drop=TRUE), function(x) rbind(NA,x,NA))
or.norep <- do.call(rbind, or.norepplusna)
nr.norep <- nrow(or.norep)

#jpeg("figure3a.jpg", width=8.5, height=5, units="in", res=800)
pdf("figure3a.pdf", width=8.5, height=5)
par(mai=c(.6,1,.4,.5))
with(or.norep, {
  plot(1:nr.norep, OR, ylim=c(.6,1.4), type='n', xaxt='n', xlab='', ylab='Odds Ratio and 95% Confidence Interval', frame.plot=TRUE)
  segments(1:nr.norep, lci, 1:nr.norep, uci, col=c("black", "black", "black", "black", "black"), lwd=3) # Need to include 5 colors because of NA rows
  #lines(c(1:nr.norep, 1:nr.norep), c(lci,uci))
  points(1:nr.norep, OR, pch=c(26,17,19,15,26), cex=1) # Need to include 5 line types because of NA rows. 26 plots blank
  xlabels <- c("Conservative", "Donations", "Vote Margin", "Revenue", "Tenure", "Unemployment")
  axis(1, c(3, 8, 13, 18, 23, 28), xlabels)
})
abline(h=1, lty=1)
abline(h=1.2, lty=2)
abline(h=0.8, lty=2)
leg.lab <- c("All Votes", "Final Bill", "Amendments")
legend(1,1.40, legend=leg.lab, pch=c(17,19,15), cex=1.2, bty="n", y.intersp=.7)
dev.off()

########################### Figure 3b #######################################

#Plot OR for Republican and Senate
or.repub <- subset(or, var=="Republican")
or.senate <- subset(or, var=="Senate")
or.repub <- rbind(or.repub, or.senate)
or.repubplusna <- lapply(split(or.repub, or.repub$var, drop=TRUE), function(x) rbind(NA,x,NA))
or.repub <- do.call(rbind, or.repubplusna)
nr.repub <- nrow(or.repub)

#jpeg("figure3b.jpg", width=6, height=6, units="in", res=800)
pdf("figure3b.pdf", width=6, height=6)
par(mai=c(.6,1,.4,.5))
with(or.repub, {
  plot(1:nr.repub, OR, ylim=c(0,46), type='n', xaxt='n', xlab='', ylab='Odds Ratio and 95% Confidence Interval', frame.plot=TRUE)
  segments(1:nr.repub, lci, 1:nr.repub, uci, col=c("black", "black", "black", "black", "black"), lwd=3)
  #lines(c(1:nr.norep, 1:nr.norep), c(lci,uci))
  points(1:nr.repub, OR, pch=c(26,17,19,15,26), cex=1.2)
  #xlabels <- c("Donations", "Conservative", "Margin", "Revenue", "Unemployment")
  axis(1, c(3, 8), labels=c("Republican", "Senate"))
})
abline(h=1, lty=1)
abline(h=10, lty=2)
abline(h=20, lty=2)
abline(h=30, lty=2)
abline(h=40, lty=2)
leg.lab <- c("All Votes", "Final Bill", "Amendments")
legend(7, 40, legend=leg.lab, pch=c(17,19,15), cex=1, bty="n", y.intersp=.8)
dev.off()


############################# Figure 4a ####################################

#jpeg("figure4a.jpg", width=6, height=6, units="in", res=800)
pdf("figure4a.pdf", width=6, height=6)
par(mai=c(1,1,.4,.5))
newdata1 <- with(amendments, data.frame(republican=1, donate0810=seq(from =0, to=166, length.out=1000), ideology=mean(ideology), unemploy=mean(unemploy), county_rev=mean(county_rev), margin=mean(margin), senate=1, tenure_total=mean(tenure_total)))
newdata1 <- cbind(newdata1, predict(model.amend, newdata=newdata1, type="link", se=TRUE))
newdata1 <- within(newdata1, {
  PredictedProb <- plogis(fit) 
  LL <- plogis(fit-(1.96*se.fit))
  UL <- plogis(fit+(1.96*se.fit))
})

newdata2 <- with(amendments, data.frame(republican=0, donate0810=seq(from =0, to=32, length.out=1000), ideology=mean(ideology), unemploy=mean(unemploy), county_rev=mean(county_rev), margin=mean(margin), senate=1, tenure_total=mean(tenure_total)))
newdata2 <- cbind(newdata2, predict(model.amend, newdata=newdata2, type="link", se=TRUE))
newdata2 <- within(newdata2, {
  PredictedProb <- plogis(fit) 
  LL <- plogis(fit-(1.96*se.fit))
  UL <- plogis(fit+(1.96*se.fit))
}) 

with(newdata1, plot(donate0810, PredictedProb, type="n", ylim=c(0,1), ylab="Probability of Pro-Industry Vote", xlab="Donations", xaxt="n"))
with(newdata1, lines(donate0810, PredictedProb, lwd=2, col="black"))
with(newdata1, lines(donate0810, UL, lty=2, col="black"))
with(newdata1, lines(donate0810, LL, lty=2, col="black"))
with(newdata2, lines(donate0810, PredictedProb, lty=5, lwd=2, col="black"))
with(newdata2, lines(donate0810, UL, lty=2, col="black"))
with(newdata2, lines(donate0810, LL, lty=2, col="black"))
xrange <- c(0,166)
ticks <- pretty(xrange)
axis(1, at=ticks, labels=c("$0", "$50,000", "$100,000", "$150,000", "$200,000"))
leg.lab <- c("Republicans", "Democrats")
legend(100, .30, legend=leg.lab, col=c("black", "black"), lty=c(1,5), lwd=2, bty="n", y.intersp=.8)
dev.off()


############################# Figure 4b ####################################

#jpeg("figure4b.jpg", width=6, height=6, units="in", res=800)
pdf("figure4b.pdf", width=6, height=6)
par(mai=c(1,1,.4,.5))
newdata1 <- with(amendments, data.frame(republican=1, donate0810=seq(from =0, to=166, length.out=1000), ideology=mean(ideology), unemploy=mean(unemploy), county_rev=mean(county_rev), margin=mean(margin), senate=0, tenure_total=mean(tenure_total)))
newdata1 <- cbind(newdata1, predict(model.amend, newdata=newdata1, type="link", se=TRUE))
newdata1 <- within(newdata1, {
  PredictedProb <- plogis(fit) 
  LL <- plogis(fit-(1.96*se.fit))
  UL <- plogis(fit+(1.96*se.fit))
})

newdata2 <- with(amendments, data.frame(republican=0, donate0810=seq(from =0, to=32, length.out=1000), ideology=mean(ideology), unemploy=mean(unemploy), county_rev=mean(county_rev), margin=mean(margin), senate=0, tenure_total=mean(tenure_total)))
newdata2 <- cbind(newdata2, predict(model.amend, newdata=newdata2, type="link", se=TRUE))
newdata2 <- within(newdata2, {
  PredictedProb <- plogis(fit) 
  LL <- plogis(fit-(1.96*se.fit))
  UL <- plogis(fit+(1.96*se.fit))
}) 

with(newdata1, plot(donate0810, PredictedProb, type="n", ylim=c(0,1), ylab="Probability of Pro-Industry Vote", xlab="Donations", xaxt="n"))
with(newdata1, lines(donate0810, PredictedProb, lwd=2, col="black"))
with(newdata1, lines(donate0810, UL, lty=2, col="black"))
with(newdata1, lines(donate0810, LL, lty=2, col="black"))
with(newdata2, lines(donate0810, PredictedProb, lty=5, lwd=2, col="black"))
with(newdata2, lines(donate0810, UL, lty=2, col="black"))
with(newdata2, lines(donate0810, LL, lty=2, col="black"))
xrange <- c(0,166)
ticks <- pretty(xrange)
axis(1, at=ticks, labels=c("$0", "$50,000", "$100,000", "$150,000", "$200,000"))
leg.lab <- c("Republicans", "Democrats")
legend(100, .30, legend=leg.lab, col=c("black", "black"), lty=c(1,5), lwd=2, bty="n", y.intersp=.8)
dev.off()


